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Abstract  of  Dissertation  Presented  to  the  Graduate  Council  of 
the  University  of  Florida  in  Partial  Fulfillment  of  the  Requirements 
for  the  degree  of  Doctor  of  Philosophy 

A FINITE  ELEMENT-DIFFERENTIAL  METHOD  FOR 
INCOMPRESSIBLE  TURBULENT  BOUNDARY  LAYER  FLOWS 

By 

Tyne-Hsien  Chang 
December  1981 

Chairman:  Chen-Chi  Hsu 

Major  Department:  Engineering  Science 

The  applicability  as  well  as  the  accuracy  and  efficiency  of  a 
finite  element-differential  method,  which  had  been  shown  to  be  a 
very  effective  method  for  laminar  flows,  is  investigated  in  detail  for 
more  complex  steady  two-dimensional  incompressible  turbulent  boundary 
layer  flow  problems.  The  closure  model  chosen  for  the  turbulent  flows 
is  a two-layer  eddy  viscosity  model.  A number  of  important  transform- 
ations have  been  carried  out  for  the  system  of  governing  equations 
before  the  application  of  the  proposed  solution  method. 

In  the  method  of  solution,  the  transformed  partial  differential 
equation  is  first  reduced  to  a system  of  first  order  nonlinear 
ordinary  differential  equations  by  a subdomain  collocation  method, 
in  which  the  unknown  function  at  a streamwise  station  is  represented 
by  a classical  spline  function.  The  resulting  initial  value  problem 
is  then  integrated  numerically  by  the  modified  Hamming's  predictor- 
corrector  method  as  well  as  by  Gear's  method  for  stiff  equations. 


vi 


The  numerical  experiments  have  been  conducted  on  the  flat  plate 
problem  which  consists  of  the  laminar,  transitional,  and  turbulent 
flow  regions  covering  the  range  of  local  Reynolds  numbers  from 
8 x 103  to  1 x 108.  The  study  shows  that  the  method  of  solution  can 
be  very  efficient  and  provide  highly  accurate  results  for  the 
turbulent  flow  problem.  In  fact  the  computed  skin-friction  coeffi- 
cient distribution  over  the  range  of  the  Reynolds  numbers  considered 
agrees  almost  perfectly  with  the  measured  data.  Moreover,  the 
results  of  numerical  experiments  have  provided  some  valuable 
information  on  the  flow  characteristics  relevant  to  numerical 
approximations . 
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CHAPTER  I 
INTRODUCTION 


In  recent  years  there  has  been  considerable  research  effort  in  the 
field  of  computational  turbulent  boundary  layer  flows  because  of  the 
need  for  more  accurate  solutions  to  flow  problems  that  arise  in  engi- 
neering applications.  The  turbulent  boundary-layer  flows  are  generally 
more  complex  than  laminar  boundary-layer  flows  in  the  sense  that  these 
flow  problems  do  have  their  own  distinct  and  important  characteristics 
which  may  require  special  attention  in  the  application  of  solution  tech- 
niques. For  instance,  the  growth  of  turbulent  boundary-layer  thickness 
occurs  at  a much  faster  but  unknown  rate,  and  the  very  thin  viscous  sub- 
layer in  the  vicinity  of  the  wall  has  an  extremely  large  velocity  gra- 
dient; consequently,  special  measures  are  needed  in  the  discretization 
of  the  flow  region.  In  general,  large  numbers  of  grid  points  with 
relatively  small  grid  sizes  must  be  specified  near  the  wall  in  order  to 
obtain  an  accurate  solution. 

The  governing  boundary  layer  equations  for  incompressible  turbulent 
flows  contain  a term,  called  the  Reynolds  stress,  which  involves  the 
time  mean  of  the  product  of  two  fluctuating  velocities.  The  Reynolds 
stress  has  been  modeled  by  means  of  the  mean  velocity,  of  the  turbulent 
energy,  of  the  product  of  turbulent  energy  with  a length  scale,  and  of 
the  Reynolds  stress  itself  as  the  dependent  variables  of  the  differential 
equations;  consequently  the  closure  models  for  turbulent  boundary-layer 
flows  can  be  classified  as  zero  equation,  one-equation,  two-equation. 
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and  stress-equation  models.  The  zero-equation  models  are  easier  to 
formulate  and  calibrate  and  therefore  have  been  used  extensively  to 
provide  predictions  for  engineering  problems.  Other  closure  models  are 
more  complex  and  their  use  for  practical  problems  is  rather  limited  at 
present.  Reynolds  (1976)  provided  an  overview  of  the  closure  models 
and  concluded  that  the  zero-equation  model  has  been  applied  to  a variety 
of  turbulent  flows  and  works  quite  will  for  most  turbulent  boundary 
layer  flows  considered  at  Stanford  conference,  a special  conference 
was  held  at  Stanford  Universtiy  in  1968  to  compare  the  viability  of 
prediction  methods  based  on  various  closure  models  for  the  partial  dif- 
ferential equations  describing  turbulent  boundary  layer  flows.  In 
addition,  the  Stanford  conference  provided  a direct  check  on  the  accur- 
acy of  the  various  prediction  methods.  A total  of  nine  methods  of 
solving  governing  partial  differential  equations  with  a closure  model 
competed.  The  results  of  the  conference  indicated  that  the  prediction 
methods  of  Mellor  and  Herring  (1968),  Cebeci  and  Smith  (1968),  and  Ng, 
Patankar  and  Spalding  (1968)  can  provide  acceptably  accurate  solution 
for  the  flow  problems  considered.  The  solution  of  Mellor  and  Herring 
with  the  effective  eddy  viscosity  model  is  obtained  completely  across 
the  boundary  layer.  They  used  the  Hartree-Womersley  semi -discrete 
method  for  solving  the  governing  equations,  where  the  ordinary  differen- 
tial equations  are  solved  across  the  layer  using  the  Runge-Kutta  inte- 
gration procedure.  Ng,  Patankar  and  Spalding  selected  a mixing  length 
model  and  used  a Crank-Nicol son  type  of  finite-difference  scheme,  which 
is  not  applied  to  the  wall.  Instead,  a wall  function  is  introduced 
that  gives  the  flow  quantities  at  the  wall.  Cebeci  and  Smith  used  an 
implicit  finite-difference  scheme  to  obtain  the  solution  across  the 
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complete  boundary  layer  with  the  eddy  viscosity  model.  Several  hundred 
grid  points  across  the  boundary  layer  had  been  used  by  the  above  men- 
tioned finite  difference  methods  for  accurate  solution. 

Recently,  the  two-layer  eddy  viscosity  closure  model  has  been  ex- 
tensively used  in  solving  practical  turbulent  flow  problems.  This 
model  consists  of  inner  and  outer  regions,  and  the  distributions  of 
eddy  viscosity  are' described  by  two  separate  empirical  expressions  in 
the  two  regions  (e.g.  Cebeci  and  Smith,  1974).  The  expression  for  eddy 
viscosity  in  the  inner  region  is  based  on  mixing  length  theory,  and  a 
modified  expression  for  the  mixing  length  is  used  to  account  for  the 
viscous  sublayer  close  to  the  wall.  The  expression  for  eddy  viscosity 
in  the  outer  region  is  based  on  a constant  eddy  viscosity  modified  by 
a transverse  intermittency  factor.  This  model  has  been  accepted  as  the 
standard  form  of  the  zero-equation  model  and  used  by  many  researchers 
to  develop  computational  methods  for  turbulent  flow  problems. 

Keller  and  Cebeci  (1972)  applied  a Levy-Lees  Transformation  to  the 
governing  equations  and  then  used  the  well  known  box  scheme  with 
Richardson  extrapolation  to  obtain  the  solution  for  turbulent  flow  over 
a flat  plate.  The  flow  region  considered  is  from  the  local  Reynolds 
number  Re  = 0.83  x 10  to  Re  = 1.88  x 10  with  constant  boundary 

X X 

layer  thickness.  The  numerical  solutions,  using  Keller's  box  scheme 
based  on  11,  21,  and  31  grid  points  across  the  boundary  layer  with  the 
designed  variable-grid  system,  were  obtained  at  17  different  stations 
in  the  streamwise  direction.  The  results  of  the  31-grid  scheme  are 
acceptably  accurate  but  the  results  of  the  other  two  schemes  are  not 
acceptable;  however,  the  accuracy  of  solutions  can  be  improved  by  the 
Richardson  extrapolation.  These  results  show  that  Keller's  box  scheme 
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with  Richardson  extrapolation  can  be  very  simple  and  accurate  for  tur- 
bulent flows. 

Blottner  (1974)  pointed  out  that  Keller's  box  scheme  requires  the 
solution  of  block- tri diagonal  equations  and  can  be  very  time  consuming 
when  extended  to  compressible  boundary  layer  flows.  Consequently,  he 
proposed  a Crank-Nicol son  type  of  difference  scheme  with  a variable 
grid  and  solved  the  same  flat  plate  problem.  The  results  obtained 
indicated  that  the  Crank-Nicol son  type  scheme  is  easier  to  formulate 
and  can  give  the  same  accuracy  as  Keller's  box  scheme  for  the  turbulent 
boundary  layer  flow.  In  addition,  the  number  of  grid  points  required 
for  the  Crank-Nicolson  method  is  less  than  that  of  Keller's  box  scheme 
under  the  same  accuracy  requirement. 

Recently  Rubin  and  Khosla  (1977)  have  extended  their  higher-order 

collocation  method  to  turbulent  flows  with  some  success.  The  higher- 

order  collocation  method  derived  from  polynomial  spline  interpolation 

and  Hermite  collocation  is  used  to  overcome  the  stiffness  problem 

caused  by  the  extremely  small  grid  size  near  the  boundary.  Again,  they 

chose  the  flow  over  a flat  plate  as  the  testing  example.  The  solutions 
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were  obtained  for  the  range  5.45  x 10  <.  Re  <.2x10.  For  the  Rey- 

X 

nolds  number  considered  here  they  only  used  11  grid  points  across  the 
boundary  layer,  with  first  element  size  Ik  = 0.1432,  to  provide  5-10% 
accuracy.  With  21  points  and  h-|  = 0.005  the  results  are  well  within 
the  experimental  scatter.  Indeed  the  method  can  also  provide  accurate 
results.  It  should  be  pointed  out,  however,  that  an  iterative  method 
can  lead  to  erroneous  solutions  if  the  truncation  error  of  the  approxi- 
mation is  large.  This  fact  is  evident  from  the  results  of  the  numerical 
experiments  given  in  Rubin  and  Khosla  (1977). 
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More  recently  a finite  element-differential  method  has  been  develop- 
ed and  tested  successfully  for  two-dimensional  laminar  boundary  layer 
flows  (Hsu,  1976  and  Hsu  and  Liakopoulos,  1981).  The  application  of 
the  method  is  very  simple  and  straightforward;  however,  one  must  carry 
out  a number  of  transformations  in  the  analysis  first.  In  this  method 
the  transformed  flow  region  is  divided  into  a number  of  strips  parallel 
to  the  boundary  and  the  unknown  functions  at  a given  streamwise  station 
are  approximated  by  classical  cubic  splines.  The  transformed  governing 
partial  differential  equations  are  then  reduced  to  a system  of  first 
order  nonlinear  ordinary  differential  equations  by  a subdomain  collo- 
cation method;  consequently,  the  resulting  initial  value  problem  is 
solved  by  a numerical  integration  method.  The  stiff  problem  associated 
with  the  reduced  equations  may  be  pronounced  if  one  used  a great  number 
of  strips  and  very  small  element  sizes  near  the  wall.  The  numerical 
results  obtained  showed  that  the  method  is  indeed  very  efficient  and 
can  provide  highly  accurate  results  for  the  entire  flow  region. 

It  seems  that  the  finite  element-differential  method  developed  is  a 
rather  promising  solution  technique  for  more  complex  flow  problems. 

The  objective  of  this  study  is,  therefore,  to  investigate  the  applicabi- 
lity as  well  as  the  accuracy  and  efficiency  of  the  finite  element-dif- 
ferential method  to  more  complex  two  dimensional  incompressible  turbulent 
boundary-layer  flows.  The  governing  boundary  layer  equations  with  the 
two-layer  eddy  viscosity  closure  model  have  been  transformed  into  the 
proper  forms  for  applying  the  solution  technique. 

The  applicability,  accuracy  and  efficiency  of  the  solution  method 
for  complex  turbulent  boundary  layer  flows  have  been  investigated  in 
detail  upon  the  flat  plate  problem.  The  boundary  layer  flow  considered 
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consists  of  the  laminar,  transitional,  and  turbulent  regions  which 

3 

covers  the  range  of  local  Reynolds  numbers  from  Re  = 8.1  x 10  to 

A 

g 

Re  = 1.1  x 10  . The  results  obtained  have  shown  that  the  finite 

A 

element-differential  method  can  be  very  efficient  and  provide  highly 
accurate  solution  for  the  turbulent  boundary  layer  flow  problems. 
Moreover,  the  result  of  numerical  experiments  has  given  valuable 
information  on  the  flow  characteristics  relevant  to  numerical  approxi- 
mation. 


CHAPTER  II 

GOVERNING  BOUNDARY-LAYER  EQUATIONS 


For  two-dimensional,  steady,  incompressible  turbulent  boundary 
layer  flows,  the  governing  equations  are  (e.g.  White,  1974) 


x = the  coordinate  along  the  boundary  in  the  direction  of  flow, 
y = the  coordinate  normal  to  the  boundary, 
u = mean  velocity  component  in  x direction, 
v = mean  velocity  component  in  y direction, 
u'  = fluctuating  velocity  component  in  x direction, 
v'  = fluctuating  velocity  component  in  y direction, 

U = velocity  just  outside  the  boundary  layer, 
y = dynamic  viscosity, 
p = fluid  density. 


The  associated  boundary  and  initial  conditions  considered  are 

u ( x , 0)  = 0 , v ( x , 0)  = 0 , u ( x , y -*  ~)  -*  U(x)  , (2-3) 


In  equation  (2-2) -pu' v'  is  called  Reynolds  stress  which  is  caused  by 
the  action  of  turbulence.  When  the  Reynolds  stress  is  zero,  equations 


(2-1) 


(2-2) 


where 


u(0,  y)  = uQ(y)  . 


(2-4) 
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(2-1)- (2-4)  are  the  governing  equations  for  laminar  boundary- layer  flows. 
Equations  (2-1 )-(2-4)  are  not  complete;  a closure  model  to  make  the 
system  closed  is  needed.  The  purpose  of  a closure  model  is  to  re- 
late the  Reynolds  stress  to  other  flow  quantities  and  thereby  obtain  a 
closed  system  of  equations. 

There  are  several  levels  of  closure  models  available  (Reynolds, 
1976).  The  zero-equation  model,  based  on  algebraic  eddy  viscosity  for- 
mulation, has  been  tested  for  a wide  range  of  turbulent  boundary  layer 
flows  with  considerable  success.  The  one-equation  model,  in  which  a 
partial  differential  equation  for  the  turbulence  energy  is  solved  in 
conjunction  with  the  partial  differential  equations  for  the  mean  motion, 
finds  use  in  practical  applications  on  occasion.  The  two-equation 
model,  obtained  by  adding  an  additional  partial  differential  equation 
for  the  turbulence  length  scale  to  the  one-equation  model,  has  not  been 
used  extensively  for  engineering  problems,  although  it  is  currently 
popular  among  researchers.  The  stress-equation  model,  which  involves 
partial  differential  equations  for  all  components  of  the  turbulent 
stress  tensor,  is  now  under  extensive  development. 

Due  to  its  simple  formulations  and  wide  applications  in  various 
engineering  problems,  the  two- layer  eddy  viscosity  model  is  used  in  the 
present  analysis.  The  development  of  this  model  is  discussed  in  some 
detail  in  appendix  A.  According  to  the  model,  the  kinematic  eddy  vis- 
cosity, v^,  is  defined  by 

_ W 


3U 

ay 


(2-5) 
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The  effective  eddy  viscosity  is  defined  by 

ve  = + v , (2-6) 

where  v is  the  kinematic  laminar  viscosity.  Consequently,  the  distri- 
bution of  the  two- layer  effective  eddy  viscosity  across  the  boundary 
layer  including  the  transitional  and  transverse  i ntermittency  factors 
can  be  written  as  (appendix  A,  equations  (A- 12)  and  ( A- 13)) 

(ve)1  = v + (0.4y[l  - exp  C-^)]H|yl  Yt  . (2-7) 

(v  )g  = v + 0.0168  U(x)S*  y > (2-8) 

where  subscripts  i and  0 represent  the  inner  and  outer  regions,  respec- 
tively. The  damping  constant  A,  transitional  i ntermittency  factor 
Yt,  transverse  intermi ttency  factor  y , and  displacement  boundary  layer 
thickness  6*,  as  discussed  in  appendix  A, are 

A = (26^  ) (1  - 11.8  , (2-9) 

T T 

* TjfsJ  ] ’ <2'10) 

tr 

Y = [1  + 5.5  (y/fi)6]"1,  (2-11) 

5*  = f (1  - -n-)dy  , (2-12) 

0 u 

where  x.  is  the  location  of  the  start  of  transition  and  6 is  the  bound- 
tr 

ary  layer  thickness.  The  friction  velocity,  u^ , and  the  transition 

Reynolds  number,  R ,are  defined  by  the  following  expressions,  respectively, 
xtr 


= 1 - exp  [ 


-1 

1200 


■1.34 


(x  - 


tr 


x,  ) 

tr 
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u 

T 


H (tjxj/p)* 
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U(x)Xtr 


(2-13) 

(2-14) 


where  t, ,(x)  is  shear  stress  at  wall, 
w 

The  matching  point  which  separates  the  inner  and  outer  regions  must 
be  determined  from  the  continuity  of  the  eddy  viscosity,  that  is. 


(ve)i  = (ve)0 


(2-15) 


Therefore,  the  complete  set  of  the  governing  boundary  layer  equations 
for  the  flow  problem  considered  is  given  by 


3u  + 3v 
3x  3y 


= 0 , 


(2-16) 


+ 


_ 3 / 3 U \ 

' ay  lveay' 


+ 4U- 

dx 


(2-17) 


1(v  ) . given  by  equation  (2-7) 

e 1 (2-18) 

( Vg ) o given  by  equation  (2-8)  , 


u(x,  0)  = 0 , v(x,  0)  = 0 , u(x,  y-**>  ) ->  U(x)  , 


(2-19) 


u(0,  y)  = Ug(y) 


(2-20) 


The  important  physical  quantity  besides  the  velocity  flow  field  is  the 
wall  shear  stress  distribution,  t (x),  which  is  defined  by 


Tw(x)  E [pve  Ty-^y=0 


(2-21) 


As  usual  the  local  skin-friction  coefficient  is  defined  by 


r _ Tw(x) 


(2-22) 


CHAPTER  III 

TRANSFORMED  BOUNDARY -LAYER  EQUATIONS 


In  solving  boundary  layer  flow  problems  governed  by  equations  (2-16)- 
(2-20),  it  is  often  advantageous  to  carry  out  a number  of  coordinate 
transformations  before  application  of  a solution  method.  First,  the 
governing  equations  can  be  made  dimensionless  by  the  following  trans- 
formati ons: 


x 

X1  L 


(3-1) 


U00 

um 


in  which  L and  Um  are  the  characteristic  length  and  velocity.  The  re- 
sulting non-dimensional  equations  are 


3 u av, 

1 + ^ = 0 , 


ax 


i 


au-j 

'l  "ax^ 


3y 


l 


+ V. 


au 


t t au, 

1 _ _o  /_  e 1 

D r\  ^ w 


) + u 


dU] 

ldx 


-yiL  3Ui 


1 ay1  Re  ay i ayn 

(Mi  = 1 + Re  {0.4  y i [ 1 - exp(— ^-)]  }|  gyL|yt 
v 1 

(^)0  = 1 + 0.0168  R5*U1(x1)Yyt  , 


(3-2) 


(3-3) 


(3-4) 


u-|(x-|,0)  = 0 , v-j  (x-| , 0)=  0 , u-j  (x-|  ,y-|-*-“) -*  U-j  (x-j ) , (3-5) 


11,(0,  y,)  = U0(y(yi))/Ujo 


(3-6) 
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where  Re  = Uoo  L/v  is  the  Reynolds  number  while  R^*  = Uoo  6*/v  is  the  Rey- 
nolds number  based  on  the  displacement  thickness  6*. 

By  use  of  another  transformation  the  transformed  equations  can  be 
made  independent  of  the  Reynolds  number  and  the  boundary  conditions  the 
same  for  all  conceivable  problems  of  the  class  considered.  Let 
x. 


x2  = J'g1  ds  ’ y2  = Reiui(x-,)y1  , 

(x-i » y-i ) 

Mx2’  = ui  (x-,)  ’ U2^x2^  = W ’ 


(3-7) 


v i ( x -j  j y-j ) Re  U-j  (x-j  )[v2(x2>  Y2)  ~ Y2  ^2  d x2  ^ " 


d In  U, 


Then  equations  ( 3-2 ) - (3-6 ) become 


3u9  3v9 

— - + — - = n 

3x2  3y2 


(3-8) 


3U?  3Up  3 

c + V 


ve  3U2 


iTr-zr)  + 0 " «4) 


2 d In  U2 


2 3x2  2 3y2  3y2  vv  3y2  v 2'  d x 


(3-9) 


= (“)  = 1 + rAo.4  y2  [1  - exp  ( — - 

v v i L ARe2U2(x2) 


(^)  = 1 + 0.0168  R5*U2(x2)YYt  , 


)]}‘ 


3U, 


Yt 

(3-10) 


u2(x2,  0)  = 0 , v2(x2’  0)  = 0 > u2(x2,  y2^°)  -*■  1 , (3-11) 


u2(0,  y2)  = uQ(y(y2))/U2(0)Uo 


(3-12) 
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It  is  known  that  the  number  of  dependent  variables  can  be  reduced 
by  one  if  the  Von  Mises  transformation  is  introduced  (e.g.  Yih,  1977). 
One  can  introduce  a stream  function  ^(x^,  y2)  to  satisfy  the  continuity 
equation  (3-8),  with 

sv'p (x2 5 ^2)  -V 2^ 

u2(x2,  y2)  = > v2(x2,  y2)  = - ^ • ( 

Then  in  the  Von  Mises  transformation  the  stream  function  4^  is  used  as 
the  independent  variable  corresponding  to  y2.  Let 


(3-14) 


u3 ( X3 , 4*3 ) - u2(x2,  y2)  , ^3^3^  U2(x2) 


Equations  (3-8)- (3-12)  become 


d In  ^^(x^) 


(3-15) 


u3%^l{0-4  y2  <*3)[1  ' exp  (' 


1 

A Re~ L x^ 7 

(3-16) 


v 


e 


= 


V 


0168  Rg*  U3(x3)wt  > 


6*  3 v 3 


U3VX3 , 0)  0 , ^3^x35  4,3-*,°^)  5 


(3-17) 


03(0,  4i2)  (y (4*3) )/iJ2 (o)Uop 


(3-18) 


Furthermore,  if  one  introduces  the  new  function 
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w 


3(x3,  1|J3)  = u3(x3,  y3) 


(3-19) 


one  obtains 


3w 


i ^3 

8x3  "3  9^3'v  9^3 


= W- 


v 3w~  d In  U^(x^) 

+ 2(1  - w3)  j J 


d x- 


(3-20) 


v , - Yo(i|io)  L 2 ^wo 

<TT>.  ' 1 + We5(0.4  y2(.,3)[1  - exp(  . - )]}  i-j^l Yt 


V 


Afyu3(x3) 


(3-21) 


fc“)  = 1 + 0-0168  R5*U3(x3)YYt  . 


VJ2^X3i  ^ 0 5 W3(x3  , ij^-**30)  "►  1 y 


(3-22) 


w3(0 , ) ^0(1^3) 


(3-23) 


For  a boundary  layer  profi1e:  whose  development  with  respect  to  y-j 
starts  with  the  linear  term,  the  independent  variable  ^ is  inconvenient 
to  use  for  a function  which  is  analytic  in  y-j  and  not  analytic  in  terms 
of  iji3  (Guderley  and  Hsu,  1972).  This  motivates  the  following  transfor- 
mati on: 


(3-24) 

W4U4,  04)  - ^3(^3’  n3)  ’ 04^4)  - U3(x3)  , 


Hence,  one  finds 

^ = wL  1 -i_(- £ ^ . 1 le  ^4,  + , d 1n  u4<x4> 

3X4  4 LoZ4  3114  V 8114  v 804^  W4'  d X4 


, (3-25) 
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v 


v 


w4(x4,  0)  = 0 , w4(x4,  n4-*»)  * 1 » 


(3-27) 


w4(0,  n4)  = w0(n4)  . 


(3-28) 


For  any  numerical  method  of  solution,  it  is  desirable  to  have  the 
transformed  boundary  layer  thickness  remain  approximately  the  same  order 
of  magnitude  through  the  entire  flow  region.  The  rate  of  growth  of  a 
turbulent  boundary  layer  is  not  known;  however,  for  laminar  boundary 
flows  the  application  of  the  Falkner-Skan  transformation  seems  to  main- 
tain a boundary  layer  thickness  of  the  same  order  throughout  the  entire 
flow  region  (Schlichti ng , 1979).  Therefore,  one  hopes  the  Falkner-Skan 
transformation  may  somewhat  decrease  the  growth  rate  of  the  turbulent 
boundary  layer.  The  Falkner-Skan  transformation  is  given  by 


where  e is  a small  constant  introduced  to  overcome  the  singularity  resulted 
from  a coordinate  transformation.  Then  equations  (3-25)- (3-28)  become 


(3-29) 


w(5,  n)  = w4(x4,  n4)  , WU)  = U4(x4)  , 
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v 


e 


(3-31) 


v 


= 


(^i)  = l + 0.0168  R6*  W (5)YYt  » 


0 


w(5,  0)  = 0 , w(s,rr*»)  -*  1 


(3-32) 


w(0,  n)  = Wg(n)  • 


(3-33) 


In  the  present  formulation,  the  dependent  variable  w(5,  n)  is  the 
square  of  the  dimensionless  velocity.  Because  of  the  homogeneous  bound- 
ary conditions,  equation  (3-11)  and  the  transformation,  equation  (3-14) 
and  (3-19),  one  has,  for  sufficiently  small  value  of  n» 


Furthermore,  the  asymptotic  behavior  of  a boundary-layer  profile  implies 
that 


For  application  of  a numerical  method,  the  boundary  condition  at 
infinity  must  be  imposed  at  a finite  but  sufficiently  large  distance  H. 

The  proper  choice  of  H is  in  general  crucual  to  the  efficiency  of  a 
given  numerical  scheme  as  well  as  to  the  accuracy  of  the  results  obtained. 
Hence,  for  a selected  H,  the  initial-boundary  value  problem  to  be 


w(5,  n)  = 0(n2)  , 


(3-34) 


which  implies  that 


(3-35) 
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numerically  solved  is  governed  by 


|f  = R,(w) 


(3-37) 


e _ 


(—)  =!+  0.25  (-%  {0.4  y2(n)[l  - 


ee 


- y^)1-  2 

exp( 1 )]} 

ARe2W(5) 


3W 


3n 


YtiT 


(-r)  = l + 0.0168  R * W(5)yyt  , 

V r\  Ci  ^ 


at  n = 0 
at  n = H 
at  e,  = 0 


w = 0 , 

w = 1 , 

w = wQ(n)  , 


(3-38) 

(3-39) 

(3-40) 

(3-41) 


where  the  nonlinear  operator  R-j  (w)  is  given  by  the  right  hand  side  of 
equation  (3-30).  The  additional  boundary  conditions,  equation  (3-35) 
and  (3-36)  are 


= 0 at  n = 0 and  at  n = H 
an 


(3-42) 


which  are  important  and  useful  in  the  spline  fit  approximation. 

Following  the  transformations  given  by  equations  (3-1),  (3-7), 
(3-14),  (3-19),  (3-24)  and  (3-29),  the  shear  stress  is  given  by 


,F>  , .(UJIsliy 

w 4Re(ee5)*  ” 3,1  0 


T 


(3-43) 


One  notes  from  equation  (3-35)  that 


32w 


/ 1 3W \ _ 

n 3n^n_Q  3nz 


n-0 


(3-44) 


Consequently,  the  local  skin-friction  coefficient,  defined  by 
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equation  (2-22),  becomes 


C 


f 


1 

2(Reee5)2 


a2w 

an7 


(3-45) 


The  system  of  governing  equations  (3-37)- (3-41)  will  be  used  to 
investigate  the  development  and  applicability  of  a finite  element  dif- 
ferential method  in  which  the  unknown  function  w(£,  n)  at  a streamwise 
station  5 is  represented  by  the  classical  cubic  spline.  Some  of  the 
important  variables  and  parameters  which  must  be  evaluated  in  the 
course  of  computation  are 


(3-46) 


y - [i  + s.se)6]-1 

nH 


(3-47) 


Yt 


= 1 


exp[- 


Re2W2(g)e2n 
l'200  Rxtr 


Ry  = Re  W(?)e 
xtr 


je^tr 

0 


ds 

WlsT  ’ 


(3-48) 

(3-49) 


- (Re  ee5)*  [J  - 2„)  d„] 
> 0 w2 


(3-50) 


CHAPTER  IV 

CLASSICAL  SPLINE  APPROXIMATION 


In  approximating  a mathematical  function,  a spline  function  will 
usually  provide  closer  and  smoother  approximations  to  the  function  and 
its  low-order  derivatives  than  a polynomial  (Greville,  1969).  The  amount 
of  computation  involved  is  not  materially  greater  than  that  required  in 
a polynomial.  In  addition,  the  spline  functions  possess  certain  impres- 
sive optimal  properties  such  as  minimum  norm,  strong  convergence  and 
best  approximation  (Ahlberg  et  al . 1967). 

A spline  function  is  a function  whose  graph  is  a composite  curve 
made  up  of  a number  of  polynomial  arcs  of  a given  degree  fitted  togeth- 
er in  such  a way  that  the  junctions  of  successive  arcs  are  as  smooth  as 
they  can  be  made  without  going  toa  si ngle  polynomial  over  the  entire 
range.  Mathematically,  the  spline  function  can  be  defined  as  follows: 
Let  0 = 0q5  n-|>  ••••  > nn>  nn+i  = H be  a strictly  increasing  sequence 

of  real  numbers  (called  nodes  of  the  spline  function).  Then  a spline 

function  of  degree  m with  the  nodes  n-j , , nn»  nn+-|  is  a 

function  s(n)  satisfying  the  following  two  conditions: 

(1)  In  each  interval  (n^ , ) (i  = 0,  1,  ...  , n)  S(n)  is  given  by 

a polynomial  of  degree  m. 

(2)  S(o)  and  its  derivatives  of  order  1,2,  m - 1 are  continuous 

over  the  region  (0,  H). 

A spline  function,  in  general,  may  be  represented  in  two  distinct, 
but  completely  equivalent  ways.  The  choice  of  the  representation 
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materially  affects  the  formulation  of  the  problem  and  its  ease  of  solu- 
tion. 

The  first  representation  is  referred  to  as  the  piecewise  poly- 
nomial representation;  it  simply  uses  different  polynomials  of  the  same 
degree  m,  one  for  each  interval.  This  formulation  can  be  written  as 

Sj(n)  = £ C.j j.(n  ~ Tij_-j)  , n e (nj_-j>  <3  ~ 2,  ••  » ^ + 

(4-1) 

S.(n)  represents  the  polynomial  of  degree  m in  the  jth  interval.  C.  . are 
J ■ J 

the  coefficients  for  the  polynomials  and  are  determined  by  applying  the 
continuity  conditions  at  internal  nodes  and  boundary  conditions  at  the 
two  extreme  nodes.  When  this  formulation  is  used,  the  continuity  con- 
ditions are  explicitly  part  of  the  problem  to  be  programmed.  This 
representation  is  easily  understood  and  the  computations  involved  are 
simple  and  straightforward. 

The  second  representation  is  formed  by  a set  of  basis  functions. 

A spline  function  represented  by  a set  of  basis  functions  is  called  a 
B-spline.  The  B-spline  may  be  written  as 


n+m+l 

S(n)  = l d.B.  (n) 
.-“l  i i ,nr  ' 


(4-2) 


where  B^  m(n)  is  the  so-called  "compact  support  basis"  and  d^  are  the 
coefficients  of  the  B-spline.  When  a B-spline  is  used  to  approximate  the 
solution  of  a differential  equation,  it  leads  to  well-behaved  matrices 
in  addition  to  having  built-in  continuity.  The  drawback  of  this  repre- 
sentation is  that  B-splines  are  not  simply  related  to  elementary  func- 
tions, and  must  be  evaluated,  using  their  definition  as  a divided 
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difference  of  the  truncated  power  function.  Further  reference  to 
B-splines  and  their  properties  is  given  elsewhere  (De  Boor,  1972). 

In  general,  cubic  spline  functions  with  m = 3 are  the  most  useful 
and  interesting  ones  because  of  their  easy  application.  The  first  repre- 
sentation of  the  cubic  spline  function  will  be  used  to  approximate  the 
transformed  velocity  profile  in  this  study. 

Cubic  spline  functions  have  been  effectively  applied  to  some  fluid 
flow  problems.  Panton  and  Sallee  (1975)  have  considered  heat  conduction 
and  laminar  boundary  layer  flow  problems  using  both  point  collocation 
and  integral  methods,  in  which  assumed  solutions  were  represented  by 
B-splines.  Rubin  and  Khosla  (1976,  1977)  have  investigated  in  detail 
the  applicability  of  cubic  splines  for  Burgers  equation,  potential 
flow  and  boundary  layer  flow  problems  using  higher-order  collocation 
methods.  Hsu  (1976)  has  studied  the  classical  cubic  spline  with  the 
subdomain  collocation  method  for  laminar  boundary  layer  flows  past  a 
submerged  body.  Furthermore,  Inoue,  Kuromaru  and  Yamaguchi  (1979) 
made  use  of  a cubic  spline  in  solving  a potential  flow  problem. 

In  the  proposed  method  of  solution  the ‘ transformed  boundary  layer 
region  is  first  divided  into  n + 1 strips  of  unequal  size  which  are 
parallel  to  the  fixed  boundary  n = 0.  The  unknown  solution  w(c,  n) 
at  a streamwise  station  E,  is  approximated  by  the  classical  cubic  spline. 
The  choice  of  the  approximation  seems  to  be  very  well  suited  for  the 
laminar  boundary  layer  flow  region  since  the  governing  equations  are  of 
second  order  in  n and  the  cubic  spline  does  satisfy  continuity  of  the 
dependent  variable,  its  slope  and  curvature  at  the  joint  nodes.  There- 
fore, it  is  expected  that  highly  accurate  results  will  be  obtained  in 
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that  region.  When  the  approximation  is  applied  to  the  transitional  and 
turbulent  flow  region,  the  slope  of  the  eddy  viscosity  which  is  not  con- 
tinuous at  the  matching  point  because  of  the  closure  model,  may  introduce 
some  difficulties  and  require  special  measures  in  approximations. 

Let  f(§,  n)  be  the  function  to  be  approximated  by  a cubic  spline 
function  in  the  interval  0 _<  n <.  H.  Assume  that  the  interval  is  divi- 
ded into  n + 1 elements  with  the  nodal  points 


0 = n0<n-|<. <n<  np+1  + H. 

Denote  the  function  f(Z,n)  and  its  first  and  second  derivatives  with 
respect  to  n evaluated  at  the  nodal  points  n-j  as 


f (5>  ’i-j ) = f i (?) » an 


-v 


U). 


a2f 


9n" 


n=n. 


= f. ' 1 


(0 


(4-3) 


Moreover,  for  a ith  element  the  local  coordinate  z and  the  element  size 

h are  defined  by 
i 


z = n 


i-1 


(4-4) 


0 <.  z <.  h . = n - - n.j_i  , for  i = 1 , . . . , n + 1 . 


Use  of  the  local  coordinate  instead  of  the  global  coordinate  simplifies 
much  of  the  formulation  and  computation.  Let  g.(5,  z)  be  the  cubic 
polynomial  approximating  f[Z,  n)  in  the  ith  element.  Then  the  continu- 
ity of  curvature  atnodesn^  -|  and  n^  indicates  that 


— 1 - [ ( h - - z) 


3Z' 


hi 


1 


fUi 


(0  + z 


f. 


(0] 


(4-5) 
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Integrating  equation  (4-5)  twice  with  respect  to  z and  applying  the 
conditions  g-U,  0)  = f^U)  and  g.(c,  h^  =»  ^.(0  one  obtains  the 
cubic  polynomial  for  the  ith  element 

[(h.  - z)3  + h 2(z  - h .)]  f\',  U)  + [z3  - h2  z]  fj'U) 

9i(5.  ^)  = — 3 1 Vhf 

(4-6) 

[(hj-  z)fi_i(5)  + z ^.(S)] 

+ S ' 

Next,  evaluating  the  slope  at  the  extreme  nodes  n0  and  nn+1  and  satis- 
fying the  continuity  of  slope  at  joint  nodes,  that  is. 


9gi fe , z) 

3Z 


39^(5.  z) 

9Z 


z=0 


(4-7) 


for  i = 1,  2,  , n,  one  obtains  the  following  system  of  (n  + 2) 

linear  algebraic  equations: 


— f 1 ' + 
3 T0 


fi' 


h.  h.,h.,,  h.  , 

1 jti  i . 1 + 1 + 1 f i i + 1 +1  .pi  i 

~ 6 fi-l  3 fi  6 l+l  h, 


,_L  * _U  f + Ii±l 

S N+i  1 Vi 


(4-8) 


^n+1  i l hn+l  , 

6 n 3 Tn+1 


n+1 


n+1 

’n+1 


+ f 


n+1 


These  equations  can  be  put  into  the  matrix  form 

[ B] { f " } = [c]{f } +{f ■ } , (4-9) 

where  both  [B]  and  [c]  are  symmetric,  tridiagonal,  (n  + 2)  by  (n  + 2) 
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constant  matrices.  The  vectors  {f11}  and  {f}  have  n + 2 elements  but 
{f1}  has  only  two  elements  fg  and  f^.  It  is  clear  that  the  inverse 
of  [B]  does  exist  and  hence  one  has 


Cf"]  = [B]'1  [c]  {f>  + [B]-1  {f  * }. 


(4-10) 


Accordingly,  one  has 


n+1 


fi'(5)=.L“ijf(5>-ei,OfO<5)+Si;ntl  f;+1(5).  (4-11) 

3 ^ 


where  the  constants  a.,  are  elements  of  [B]*^[c]  and  3.  . are  elements 

1 J 1 5 J 

of  the  inverse  of  [B] . 

In  accordance  with  equation  (3-42)  if  one  assumes  that  the  slopes 
at  the  extreme  nodes  ng  and  nn+-j  are  zero,  one  obtains 
n+1 

f\'U)  = 1 a,,f.(c)  , i = 0,  1,  ....  , n,  n + 1.  (4-12) 

1 j=0  1J  J 


This  implies  that  the  second  derivative  at  each  node  depends  linearly 
upon  all  the  nodal  functions  f.(c);  however,  f'. '(?)  is  generally  influ- 
enced  by  the  nodal  functions  f.(§)  in  the  neighborhood  of  the  ith  node. 

J 

For  example,  when  the  element  size  distribution  is  the  following: 


5 , n 

= 3 , h. 

= /0.2,  0. 

3,  0.5,  1 

.5/  , 

determined  and 

written  in 

a matrix 

form 

as 

-101 .76 

115.99 

-15.86 

1.73 

-0.11 

53.51 

-81.98 

31.71 

-3.46 

0.22 

-10.54 

29.28 

-28.47 

10.38 

-0.65 

1.62 

-4.50 

8.07 

-7.14 

1.95 

-0.81 

2.25 

-4.04 

4.90 

-2.31 

(4-13) 
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which  clearly  indicates  the  weight  of  f.(c)  on  the  curvature  fi'U). 
Substituting  equation  (4-12)  into  equation  (4-6)  gives 

n+1 

g,(S,  z)  = l a,,(z)  f.(?)  , (4-14) 

1 j=0  J J 


where  a.,  (z)  is  a known  polynomial  of  degree  three  in  z and  is  given 
■ J 

by 

c<hi  - z>3  + h?<z  - ai-i„i + tz3  - hi2zHi 


au(z)  = 


6 h. 


hi  " z z 

+ hi  6i-l  ,j  + hi  6i  ,j  ’ 


(4-15) 


Here  <5  . .is  the  Kronecker  delta.  Therefore,  the  classical  spline 
approximating  the  function  f(?,  n)  is 

n+1 

s(?,  n)  = l 5.  g,(e,  z)  , (4-16) 

i=l  1 1 


where  6-  = 1 when  n^_-|  s n s nn- , otherwise  <$..  =.0. 


CHAPTER  V 

METHOD  OF  SOLUTION 


Natural  variational  principles  do  not  exist  for  boundary  layer 
flow  problems  and  hence  formulation  of  the  finite  element  method  should 
be  based  on  the  governing  equations.  This  approach  is  categorized  as 
the  weighted  residual  method.  The  weighted  residual  method  is  a mathe- 
matical means  of  finding  approximate  solutions  to  differential  equa- 
tions. One  of  its  most  attractive  features  is  that  this  method  is 
applicable  to  nonlinear  and  non-self-adjoint  problems.  For  a given 
differential  equation  and  associated  boundary  and  initial  conditions, 
the  general  approach  of  the  weighted  residual  method  is  to  assume  a 
trial  solution  whose  functional  dependence  is  chosen,  but  includes  un- 
determined functions.  The  latter  are  found  by  requiring  that  the  trial 
solution  satisfy  the  differential  equations  in  some  specified  approxi- 
mate sense. 

In  the  proposed  numerical  method  for  solution  of  equations  (3-37)- 
(3-41),  the  boundary  layer  flow  problem,  a subdomain  collocation  method, 
one  of  the  weighted  residual  methods,  is  employed  to  reduce  the  govern- 
ing partial  differential  equations  into  a system  of  first  order,  non- 
linear, ordinary  differential  equations.  The  reduced  initial  value 
problem  is  then  numerically  integrated  to  provide  the  solution. 

For  the  spline  approximation  of  w (5,  n)  at  a given  streamwise  sta- 
tion £,  the  interval  0 £ n <.  H is  divided  into  n + 1 elements  with  nodes 
at  0 = rig  < n-j  < c nn  < nn+-|  = H and  element  size  h^  = - n^_-| 
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for  i=l,  , n+1 . If  one  denotes  the  unknown  function  w(5,  n),  and 

its  first  and  second  derivatives  at  the  ith  node  as 


w(5,  ,,)■«,(«)  , fn 


• 92w 

= wl  (5>  ■ 3^ 


= wj'  (5)  , (5-1) 


n=n^ 


n=n. 


and  also  lets  w.(5>  z)  be  the  cubic  polynomial  approximating  w(5,  n)  in 
the  ith  element,  then  the  assumed  solution  for  the  weighted  residual 
is  the  sum  of  the  cubics  and  can  be  written  as 
n+l  _ 

w(S,  n)  = l «i  w.(g,  z)  , (5-2) 

i=l 


where 


1,  when  n^i  < n < ni 
0,  otherwise 


z = n 


M-l 


As  a result  of  spline  approximation,  equations  (4-5)- (4-12) , the  assumed 
solution  can  be  written  as 

n+l  n+l 

w(e,  n)  = l ( j 6-  a . . ( z ) ) w .(5)  , (5-3) 

j=o  i=i  1 J 


in  which  a..(z)  is  given  by  equation  (4-14).  The  quantities  inside  the 
^ J 

parentheses  are  called  shape  functions,  which  are  known  polynomials  of 
degree  three  in  z. 

The  substitution  of  the  assumed  solution,  equation  (5-3),  into  the 
governing  equation  (3-37)  and  application  of  the  weighted  residual 
method  gives 


(5-4) 
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where  y-(n)  are  the  weighting  functions  and  the  nonlinear  operator 
>1 

R-|  (w)  is  given  by  the  right  hand  side  of  equation  (3-30).  Because  of 

the  boundary  conditions  given  by  equation  (3-32),  there  are  only  n 

values  of  w.(0  for  j = 1,  , n in  equation  (5-3)  to  be  determined. 

J 

Thus  one  needs  n equations  for  n unknowns;  consequently,  it  requires  a 
choice  of  n weighting  functions  to  be  used  in  equation  (5-4).  The  error 
made  by  the  spline  approximation  at  each  interior  nodal  point  is  mostly 
propagated  to  the  neighboring  intervals  (Hildebrand,  1974).  This  sug- 
gests that  the  weighting  functions  y.(n)  may  be  chosen  so  that  the  resi- 

J 

duals  will  be  distributed  to  both  sides  of  the  nodal  point  in  some  sense. 
For  simplicity  the  weighting  functions  employed  in  this  study  are  boxcar 
functions  of  the  form 

f 1 , when  n,  ->  ^ n ^ n,-+1 

r.(n)  = J_1  J 1 (5-5) 

J I 0,  otherwise  . 


This  choice  implies  that  the  subdomain  collocation  method  is  employed 
in  the  approximate  procedure  which  forces  the  residual  to  be  uniformly 
distributed  over  the  two  neighboring  elements  for  each  interior  nodal 
point.  Equation  (5-4)  incorporated  with  equations  (5-2)  and  (5-5), 
therefore,  gives  the  following  system  of  n equations 


h,  aw.  h.  , 3w.+1  h.  h.  , 

/’  P dz  + /1+1  -jP  dz  = ' R,{w,)dz  + f1+1  R,(w  )dz,  (5-6) 
0 95  0 35  0 6 1 1 ' 


for  i = 1,  2,  ...  , n.  The  substitution  of  the  explicit  expression  for 
the  cubic  polynomial  w .(5,  z)  in  equation  (4-14)  then  gives  a system  of 
n first  order  nonlinear  ordinary  differential  equations  in  w-|(0,  , 
, wn(?)>  which  can  be  expressed  in  matrix  form  as 
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[Q]  fjf>  - (r/?.  w)>  (5-7) 

or 

t||>-  [Q]"1 {?,({,«)>  (5-8) 

where  [Q]  is  a nonsingular  n by  n constant  matrix  with  its  elements  given 
by 


qi,j  = "24^-hi  (ai , j+1  + ai+l,j+l)  + hi+l (ai+l , j+1  + ai+2,j+l):i 


h. 

+ 6 • . -?r-  + 6 . . 

1-1  ,J  2 i,J 


<hi  + hui> 


+ 5 


hi+l 


i+1 » j 2 


(5-9) 


The  elements  of  the  n dimensional  vector  {w}  are  w-|(5),  , •••  , 

wn(?) , and  those  of  {y-j  (C »w) } are  the  integrals  of  nonlinear  functions 
of  the  elements  of  {w}  given  by  the  right  hand  side  of  equation  (3-30), 
that  is, 


w 


R-)  (w)  - ^ 


1 a_,_e  3W< 

ry2"  3ri  v 3n 


1 aw  ve v -i  n aw  %d  In  W (?) 

~ — — ) J + 4 ^ + " W1 h7 


n an  v 


d5 


(5-10) 


With  the  given  initial  condition,  equation  (3-41),  which  is  discretized 
at  the  interior  nodes,  the  system  of  n equations,  (5-8),  can  be 
integrated  numerically  to  provide  the  solution. 


CHAPTER  VI 

COMPUTATIONAL  METHODOLOGY 


The  computational  work  involved  in  the  proposed  method  appears 
very  complicated  and  laborious.  However,  the  computations  can  be  sim- 
plified and  reduced  to  a minimum  by  developing  the  numerical  algorithm 
carefully. 

For  a given  boundary  layer  flow  problem,  the  constant  matrices  [B], 
[C]  and  [Q]  in  equations  (4-9)  and  (5-7)  are  determined  after  the  flow 
region  has  been  discretized  into  n + 1 elements  in  the  n direction.  The 
question  of  finding  the  inverse  of  the  matrices  [B]  and  [Q]  is  of  pri- 
mary importance.  There  are  many  methods  whereby  the  inverse  of  a matrix 
can  be  determined.  In  this  study  the  Gauss-Jordan  elimination  method 
has  been  used,  since  the  difference  in  computer  time  between  the  Gauss- 
Jordan  elimination  method  and  the  Gauss  elimination  method  is  insigni- 
ficant in  finding  the  inverse  of  matrices  such  as  [B]  and  [Q].  For  the 
mathematical  basis  of  the  Gauss-Jordan  elimination  and  a FORTRAN  program 
see  Hornbeck  (1 975) . 

The  main  difference  between  the  proposed  method  and  the  finite  dif- 
ference method  is  that  integrations  are  required  using  the  finite  element 
approach.  The  right  hand  side  of  equation  (5-7)  must  be  evaluated  at 
each  integration  step.  Unfortunately,  an  explicit  expression  for  ele- 
ments of{r-|(s,  w ) } cannot  be  obtained  due  to  the  fact  that  the  opera- 
tor R-j(w)  involves  the  square  root  of  the  dependent  variable.  Hence, 
an  approximation  technique  must  be  employed  for  evaluating  the  integrals 
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on  the  right  hand  side  of  equation  (5-6).  Two  different  approximation 
techniques  had  been  investigated  for  the  laminar  boundary  layer  flow 
problem  (i.e.  ^ = 1)  (Hsu,  1976).  The  first  approximation  utilizes  a 
Taylor  series  expansion.  The  manipulation  involved  in  this  approxima- 
tion is  rather  tedious  and  cumbersome;  hence,  this  approach  is  not 
suitable  for  the  complex  turbulent  boundary  layer  flows.  An  alternative 
approach  considered  was  to  evaluate  those  integrals  numerically  using  a 
Gauss-Legendre  quadrature  formula, whi ch  is  generally  capable  of  supplying 
comparable  accuracy  with  fewer  ordinates  compared  with  numerical  inte- 
gration formulas  such  as  the  Simpson's  rule  and  the  trapezoidal  rule. 

This  latter  approach  is  employed  in  the  present  study. 

The  Gauss-Legendre  quadrature  formula  is  given  by 

= I Af>  F(4m>)  , (6-1) 

m=l 

(Ml  (M) 

where  the  set  of  weighting  coefficients  and  roots  and  ^ , for  a 

given  M are  obtained  in  such  a way  that  the  quadrature  formula  is  exact 
whenever  the  integrand  F(?)  is  a polynomial  of  degree  less  than  or 
equal  to  2M  - 1.  For  reference,  see  Krylov  (1962).  For  actual  evalua- 
tions of  the  integrals  given  in  equation  (5-6),  one  notes  that  the  sim- 
ple transformation  of  the  variable 


hi 

z = — ~ (1  + ?) 


(6-2) 


gives  the  relation 


h.  _ h.  1 h,  M 

J 1 F(z)  dz  = -4-  j F(c)  d? 

0 c -1  " m=l 


4 1 A.W  F(,W) 

2 m x m ' 


(6-3) 
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For  turbulent  boundary-layer  flow  problems  the  values  of  M = 6 and 
M = 4 have  been  tested.  The  results  obtained  using  M = 6 and  M = 4 are 
almost  the  same.  This  implies  that  M = 4 can  give  acceptable  accuracy. 

When  one  comes  to  evaluate  the  functional  values  F(^'^)  another 
difficulty  arises  because  F is  not  an  explicit  function  of  z,  due  to 
the  equation  (3-46).  Consequently,  a simple  tabulation  and  interpola- 
tion technique  must  be  employed  in  determining  the  value  of  y^  f°r  a 
given  n at  each  streamwise  station  g.  The  interpolation  process  is 
based  on  the  polynomial  of  an  assigned  degree  constructed  by  points 
which  are  close  to  the  value  for  which  interpolation  is  desired.  This 
interpolation  is  carried  out  by  a general  purpose  function  subprogram 
ATKN  (X,  Y,  N,  k,  XI),  which  is  listed  in  appendix  B. 

The  reduced  initial  value  problem,  equation  (5-8),  can  be  solved 
by  a numerical  integration  method.  Hamming's  4th  order  predictor-cor- 
rector method  (Ralston,  1960)  and  Gear's  stiff  method  (Gear,  1971b) 
are  employed  to  solve  equation  (5-8)  in  this  study. 

The  solution  obtained  by  Hamming's  4th  order  predictor-corrector 
method  shows  that  the  integration  step  size  is  quite  small  so  that  the 
efficiency  of  the  solution  method  is  deteriorated.  This  motivates  the 
investigation  on  the  degree  of  stiffness  associated  with  the  reduced 
initial  value  problem.  Because  of  the  nonlinearity  of  equation  (5-8), 
one  has  to  examine  the  results  obtained  using  the  Hamming  method  to 
determine  the  degree  of  stiffness  of  the  equations.  The  reduced  initial 
value  problem  is  indeed  very  stiff.  The  typical  step  size  required  is 
extremely  small  for  the  solution  behavior  obtained.  Accordingly, 

Gear's  stiff  method  which  has  been  shown  to  be  very  efficient  for  stiff 
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equations  (Hindmarsh,  1974)  is  utilized  to  solve  the  reduced  initial 
value  problem. 

The  University  of  Florida  Amdahl  470  V/6-II  and  IBM  3033N  computers 
have  been  used  to  perform  the  numerical  experiments  in  double  precision 
computations.  Cases  I to  IV  have  been  performed  by  the  Amdahl  470  V/6- 
II.  Cases  V to  XII  have  been  performed  by  the  IBM  3033N. 


CHAPTER  VII 

FINDINGS  OF  PRELIMINARY  NUMERICAL  EXPERIMENTS 

To  investigate  the  applicability  as  well  as  the  accuracy  and  ef- 
ficiency of  the  solution  method  for  turbulent  boundary-layer  flows, 
the  turbulent  boundary  layer  flow  over  a flat  plate  has  been  considered 
for  numerical  experiments  in  this  study.  This  flow  problem  is  a severe 
test  for  a numerical  method,  since  the  last  term  of  the  equation  (3-30) 
is  zero  and  consequently,  the  only  forcing  function  which  is  related  to 
the  eddy  viscosity  for  changing  the  solution  profile  in  the  streamwise 
direction  is  implicitly  dependent  upon  the  unknown  solution.  For  such 
a parabolic  problem,  an  accurate  initial  profile  is  very  essential  for 
accurate  results.  Any  error  made  in  the  computed  effective  eddy  viscos- 
ity may  require  a large  number  of  integration  steps  to  correct;  moreover, 
an  inaccurate  computation  for  the  effective  eddy  viscosity  distribution 
across  the  boundary  layer  at  a streamwise  station  can  be  detrimental  to 
the  solution.  Therefore,  to  numerically  solve  the  governing  equations 
for  turbulent  flows  over  a flat  plate  with  a two-layer  eddy  viscosity 
model  is  a great  challenge. 

In  the  development  of  the  solution  method  a problem  with  a unit 
Reynolds  number  of  1 x 10^(Uoo  = 160  ft/sec,  v = 1.6  x 10  ^ ft^/sec, 

L = 1 ft)  has  been  chosen  for  the  numerical  test.  The  transitional 
point  from  laminar  to  turbulent  flows  is  specified  at  the  local  Reynolds 
number  Rex  = ~ = 5.45  x 10^  which  had  been  used  by  Rubin  and  Khosla  (1977) 
and  corresponds  to  the  transformed  variable  = £ = 10.906.  The 
transitional  intermittency  factor  for  flow  over  a flat  plate  (W(c)  = 1) 
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can  be  obtained  from  equation 
uation  (3-48)  is  equal  to  10 
is  given  by 

Yt  = 1 - exp[- 


(3-48)  under  the  assumption  that  e in  eq- 
. The  transitional  intermittency  factor 


(e 


- e^V 


1200(e  tr) 


1.34 


(7-1) 


For  application  of  the  method,  a boundary  condition  at  infinity  must 
be  imposed  at  a finite  but  sufficiently  large  distance  H from  the  fixed 
boundary.  The  interval  0 £ n £ H is  then  divided  into  a limited  number 
of  elements,  N,  of  different  sizes,  h^  for  i = 1,..,N,  for  spline 
approximations.  How  to  discretize  the  region  considered  for  a finite 
element  method  depends  completely  on  how  much  is  known  about  the  problem. 
In  turbulent  boundary  layer  flows  it  is  known  that  very  steep  velocity 
gradients  occur  near  the  wall  and  approximately  zero  velocity  gradients 
exist  near  the  edge  of  the  outer  boundary  layer.  From  the  efficiency 
view  point  of  the  method,  a variable  grid  system  is  preferred.  Thus, 
small  element  sizes  must  be  used  in  the  vicinity  of  the  wall  in  order  to 
obtain  an  accurate  estimation  of  the  wall  shear  stress,  and  relatively 
large  element  sizes  can  be  used  near  the  edge  of  the  boundary.  Since 
the  nature  of  velocity  profiles  in  turbulent  boundary  layers  is  not  well 
understood,  the  importance  and  difficulty  in  choosing  the  proper  or  op- 
timal discretization  of  the  flow  region  is  greatly  enhanced  in  this  study. 

Since  the  initial  profile,  in  general, is  not  available  in  turbulent 
regions,  the  integration  must  be  started  at  a streamwise  station  in  the 
laminar  region  where  the  initial  profile  can  be  obtained.  An  approxi- 
mate profile  to  the  Blasius  equation  is  used  as  the  initial  profile 
in  this  study.  This  transformed  approximate  profile  is 
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w0(n)  = [2(g)2  - 2(1)3  + (ft)4]2  . (7-2) 

For  the  first  numerical  experiment,  the  following  finite  element 
model  was  selected  based  on  experience  from  computation  of  laminar 
boundary  layer  flows: 

case  I:  H = 7.0,  N = 16,  h.  = /0.1,  0.2,  0.3,  0.4,  12*0.5/. 

The  value  of  H = 7.0  is  believed  to  be  large  enough  for  the  range  of  the 
local  Reynolds  number  to  be  covered.  The  integration  of  the  reduced 
initial  value  problem,  equation  (5-8),  is  started  at  a streamwise 
station.  Re  = 4.8  x 104,  in  the  laminar  flow  region  with  the  assumed 

X 

profile  given  by  equation  (7-2)  as  the  initial  profile.  The  average 

integration  step  required  for  an  average  accuracy  of  10  ^ in  the  modi- 

-4 

fied  Hamming's  4th  order  predictor-corrector  method  is  A£  = 4 x 10  . 

The  computed  local  skin-friction  coefficient,  C^, is  shown  in  Figure  1 
as  a function  of  Re  .The  CPU  (Central  Processing  Unit)  time  spent  on 

X 

this  case  is  570  seconds.  The  measured  skin-friction  coefficients  of 
Dhawan  (1952),  Smith  and  Walker  (1959),  and  Winter  and  Gaudet  (1970)  are 
also  given  in  Figure  1 for  comparison.  It  is  clear  that  the  computed 
C ^ is  in  excellent  agreement  with  the  experimental  data  up  to 

5 

Re  = 7.0  x 10  . This  implies  that  the  solution  method  is  applicable 

X 

to  complex  turbulent  flows;  however,  the  accuracy  of  the  result  starts 
to  deteriorate  further  downstream.  A close  examination  of  the  computed 
results  has  shown  that  the  accuracy  of  the  computed  effective  eddy 
viscosity  distribution  across  the  boundary  layer  starts  to  degenerate 
by  exhibiting  oscillatory  characteristics  in  the  inner  region  somewhere 
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5 

after  Re  = 3.0  x 10  ; but  a short  distance  downstream  from  Re  = 

X X 

5 

3.0  x 10  the  wiggle  occurring  in  the  effective  eddy  viscosity  has 
negligible  effect  upon  the  curvature  of  the  solution  profile  at  the 
boundary,  Wq'(s),  which  is  linearly  proportional  to  in  this  case. 
Figure  2 presents  the  distribution  of  the  effective  eddy  viscosity, 

— , computed  at  a number  of  streamwise  stations.  The  deterioration  of 
the  local  skin-friction  coefficient  may  be  caused  by  the  erroneous 
distribution  of  the  effective  eddy  viscosity  which  may  be  resulting 
from  the  distorted  solution  profiles  obtained  in  the  problem.  It  is 
hard  to  recognize  if  this  is  so  at  this  moment. 

To  gain  further  insight  onto  the  accuracy  and  efficiency  of  the 
solution  technique,  the  following  refined  finite  element  model  was 
considered: 

case  II:  H = 7.0,  N = 16,  hi  = /3*0.1,  2*0.2,  3*0.3,  0.4,  4*0.5,  3*1.0/. 

The  numerical  integration  of  the  reduced  initial  value  problem  is  started 

5 

at  Re  = 6 x 10  with  the  initial  profile  obtained  by  interpolating  the 

results  of  case  I.  The  integration  step  of  Hamming's  method  with  the 

-4 

same  accuracy  as  case  I is  about  A £=  2 x 10  which  is  only  half  of  the 
step  size  allowed  in  case  I.  This  reflects  the  increase  in  the  stiff- 
ness of  the  reduced  system  of  ordinary  differential  equations  due  to 
small  element  sizes  used  near  the  wall;  consequently,  the  efficiency 
of  the  method  is  reduced.  The  computation  for  this  case  was  stopped  at 
Rex  = 10^.  The  computed  C^.  is  also  shown  in  Figure  1 and  is  quite 
accurate.  The  computed  effective  eddy  viscosity  distributions  also 
exhibit  some  oscillatory  behavior  in  the  inner  region;  however,  the 
amplitude  of  the  oscillation  is  smaller  than  that  in  case  I at  the  same 
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local  Reynolds  number.  This  means  that  the  small  element  sizes  used 
near  the  wall  may  have  some  influence  on  the  computed  viscosity  dis- 
tribution. The  CPU  time  used  in  this  case  is  200  seconds. 

The  numerical  results  obtained  from  cases  I and  II  seem  to  indicate 
that  the  given  value  of  H may  be  too  large  for  the  flow  region  consid- 
ered. In  order  to  investigate  the  importance  of  choosing  a proper  value 

of  H we  have  also  considered  the  following  model: 

case  III:  H = 5.0,  N = 10,  h.  = /2*0.1,  0.2,  0.3,  2*0.4,  0.5,  3*1.0/. 

The  value  of  H = 5.0  corresponds  approximately  to  the  boundary  layer 

thickness  at  Re  = 2 x 10  . Again  the  initial  value  problem  is  started 

A 

at  Re  = 6 x 10^  using  Hamming's  method  with  the  initial  profile  obtained 

X 

from  the  result  of  case  I.  The  computation  for  this  case  is  terminated 

at  Rex  = 2 x 10^  and  the  CPU  time  required  is  about  200  seconds.  The 

computed  C^  shown  in  Figure  1 is  apparently  very  acceptable.  Figure  3 

presents  the  distribution  of  — computed  at  a number  of  local  Reynolds 

numbers.  Again  the  oscillatory  behavior  starts  to  show  up  after 

Re  = 1.0  x 10^.  The  efficiency  of  the  method  can  be  improved  by  a small 
X * 

value  of  H;  consequently,  a scheme  which  starts  with  small  H at  relatively 
small  boundary  layer  thickness  and  increases  automatically  with  increas- 
ing boundary  layer  thickness  would  be  a most  effective  way  for  computing 
the  turbulent  flows. 

The  comparison  between  the  computed  local  skin-friction  coefficient 
for  case  I and  case  III  shown  in  Figure  1 seems  to  imply  that  the 
deterioration  of  the  local  skin-friction  coefficient  in  case  I might  be 
due  to  the  accumulation  of  errors,  which  resulted  from  the  extremely 
large  number  of  integration  steps  used.  The  following  finite  element 
model,  which  is  the  same  as  case  III,  is  adopted  to  find  out 
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whether  the  accumulation  error  is  pronounced  or  not: 

case  IV:  H = 5.0,  N = 10,  h^.  = /2*0.1,  0.2,  0.3,  2*0.4,  0.5,  3*1.0/. 

4 

The  computation  for  this  case  is  carried  out  in  the  region  of  4.8  x 10 
< Re  < 2 x 106.  Equation  (7-2)  is  used  as  the  intital  profile.  The 

X 

computed  is  shown  in  Figure  4.  The  local  skin-friction  coefficient 
obtained  in  case  III  is  also  shown  in  Figure  4 for  comparison.  The 
difference  between  cases  IV  and  III  in  the  local  skin-friction  coeffi- 

r C 

cient  within  the  range  of  6 x 10  <.  Re  <2x10  is  insignificant. 

This  implies  that  the  accumulation  error  is  not  pronounced.  Accordingly, 
the  deterioration  of  the  local  skin-friction  coefficient  and  the  wiggle 
occurring  in  the  computed  effective  eddy  viscosity  may  result  from 
other  sources.  The  CPU  time  used  in  this  case  is  about  400  seconds. 

The  numerical  results  obtained  from  case  IV  show  that  the  solution 
method  is  indeed  applicable  to  complex  turbulent  flows  and  the  computed 
local  skin-friction  coefficient  is  in  excellent  agreement  with  the 

g 

experimental  data  up  to  a local  Reynolds  number  of  Rex  = 2 x 10  . How- 
ever, the  solution  method  is  not  efficient  because  of  the  small  inte- 
gration step  required.  This  motivated  the  investigation  on  integration 
methods  to  overcome  the  numerical  instability.  The  numerical  results 
obtained  by  Hamming's  method  indicated  that  the  reduced  equations  are 
stiff  due  to  the  fact  that  the  average  integration  step  size  required  is 
much  smaller  than  that  needed  by  the  solution  behavior.  Gear's  method 
for  stiff  equations  is  chosen  to  efficiently  solve  the  reduced  initial 
value  problem.  The  DIFSUB  subroutine  subprogram  (Gear,  1971a)  of 
Gear's  method  has  been  modified  and  incorporated  into  the  main  pro- 
gram to  solve  the  reduced  initial  problem.  The  accuracy  of  the  method  is 
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assigned  as  10"4.  Moreover,  since  the  boundary  layer  thickness  of 
turbulent  flows  increases  with  5 and  its  growth  rate  is  hard  to  foresee, 
it  is  generally  impractical  to  select  a large  constant  value  of  H for 
the  flow  region  considered.  Therefore,  the  program  is  also  implemented 
to  automatically  adjust  the  boundary  layer  thickness,  H,  basedjon  an 
empirical  formula.  In  this  study  the  boundary  layer  thickness  is 
assumed  to  be  8.5  times  the  displacement  boundary  layer  thickness.  The 
number  of  elements  is  increased  in  the  program  to  account  for  the 
growth  of  the  boundary  layer.  The  maximum  allowable  integration  step  is 
limited  to  10"^  in  DIFSUB  subroutine. 

The  following  finite  element  model  was  used  to  start  the  intital 
value  problem  for  demonstrating  the  efficiency  of  Gear's  method: 

case  V':  H = 4.0,  N = 8,  h.  = /2*0.1,  0.3,  3*0.5,  2*0.1/  . 

The  computation  is  started  right  after  the  transitional  point  at 
4 

Re  = 6.0  x 10  with  the  initial  condition  given  by  equation  (7-2)  and 

A 

is  carried  out  up  to  Re  = 2.24  x 10^.  The  final  boundary  layer  thick- 

A 

ness  and  the  number  of  elements  are  H = 6 and  N = 10,  respectively. 

The  CPU  time  spent  on  this  case  is  only  22  seconds,  which  is  approximate- 
ly 20  times  less  than  that  required  for  case  IV.  It  is  clear  that  the 
efficiency  of  the  solution  method  has  been  drastically  improved  by  using 
Gear's  method.  Moreover,  the  computed  results  showed  that  the  deterior- 
ation of  the  effective  eddy  viscosity  starts  somewhere  after  Re  = 1.1  x 106 

X 

and  the  wiggle  in  the  effective  eddy  viscosity  occurs  in  the  neighborhood 
of  n = 1.0  in  this  case;  on  the  contrary,  the  deterioration  and  wiggle 
of  the  effective  eddy  viscosity  occurs  after  Re  = 1.5  x 10  and  n = 0.7, 

X 

respectively,  in  case  IV.  When  comparing  the  element  size  distribution 
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in  case  IV  and  case  V,  one  may  tell  that  the  deterioration  of  the 
effective  eddy  viscosity  takes  place  at  different  local  Reynolds  number 
and  the  wiggle  in  the  computed  effective  eddy  viscosity  occurring  at 
different  locations  in  the  n direction  may  be  due  to  the  erroneous 
solution  profiles  which  result  from  truncation  errors.  This  fact  seems 
to  imply  that  the  effective  eddy  viscosity  can  be  very  sensitive  to  the 
solution  profiles. 

Since  the  stiffness  problem  induced  by  using  the  small  element 
sizes  near  the  wall  has  been  overcome  by  using  Gear's  method  for  stiff 
equations,  how  the  effective  eddy  viscosity  is  influenced  by  very  small 
elements  near  the  wall  can  be  investigated  numerically.  The  finite  ele- 
ment model  case  VI  was  arranged  for  this  purpose: 

case  VI:  H = 4.0,  N = 9,  h.  = /2*0.05,  0.1,  0.3,  3*0.5,  2*1.0/. 

The  computation  covers  the  region  of  6.0  x 10^  <.  Re  <.2.57  x 106. 

X 

Deterioration  of  the  effective  eddy  viscosity  takes  place  further  down- 

g 

stream  after  Re  = 2.1  x 10,  and  the  wiggle  which  occurs  in  the  com- 

A 

puted  effective  eddy  viscosity  is  still  in  the  neighborhood  of  n = 1.0. 
These  results  show  that  the  very  small  element  sizes  used  near  the  wall 
can  defer  the  occurrence  of  the  deterioration  of  the  effective  eddy 
viscosity. 

The  next  problem  one  needs  to  investigate  is  whether  or  not  the 
element  sizes  not  in  the  immediate  vicinity  of  the  wall  have  any 
influence  on  the  calculation  of  the  effective  eddy  viscosity.  The 
following  finite  element  model  was  designed  for  this  purpose: 
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case  VII:  H = 5.0,  N = 11,  h..  = /2*0.05,  0.1,  0.2,  2*0.3,  2*0.5,  3*1.0/. 

The  element  sizes  are  slightly  changed  within  the  range  of  n = 0.0  to 
n = 1.0  compared  with  case  VI.  The  computation  for  the  reduced  initial 
value  problem  is  started  at  Rex  = 2.0  x 10^  and  proceeds  to  Rex  = 1.47 
x 107.  The  deterioration  of  the  effective  eddy  viscosity  still  occurs; 
however,  the  wiggle  in  the  effective  eddy  viscosity  is  shifted  farther 
from  the  wall  in  the  neighborhood  of  n = 1.50  instead  of  n = 1.0  in 
case  VI. 

Since  in  cases  VI  and  VII  we  found  that  small  elements  near  the  wall 
can  defer  the  occurrence  of  the  deterioration  of  the  effective  eddy 
viscosity,  and  small  elements  used  not  in  the  immediate  vicinity  of  the 
wall  can  shift  the  location  of  the  wiggle  farther  away  from  the  wall, 
it  is  believed  that  a slight  variation  of  the  solution  characteristics 
can  very  significantly  affect  an  accurate  calculation  of  the  effective 
eddy  viscosity.  In  order  to  gain  further  insight  the  variation  of  the 
solution  w (?,  n),  its  first  derivative  at  n-j  and  second  derivative  at 
rig  obtained  from  case  VI  are  shown  as  functions  of  Rex  in  Figures  5 and 
6.  Upon  close  examination  of  the  results  in  these  figures,  it  is  evi- 
dent that  the  first  element  size  smaller  than  0.05  should  be  used  if  the 
computation  proceeds  to  local  Reynolds  number  greater  than  10^. 

Most  of  the  finite  difference  methods  available  for  turbulent  flows 
have  been  tested  only  for  the  region  in  which  the  local  Reynolds  number 

r 

is  less  than  2 x 10  . However,  in  this  study  we  are  especially  interes- 
ted in  applying  the  finite  element-differential  method  to  sufficiently 
high  local  Reynolds  number  in  order  to  find  some  information  on  the  flow 


character!'  sties. 
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The  Gauss-Legendre  quadrature  formula  worked  very  well  in  the 

v 

laminar  region  (i.e.  — = 1)  and  its  applicability  in  the  turbulent 
region  is  expected  to  be  good  as  well.  However,  due  to  the  two-layer 
eddy  viscosity  model,  in  which  the  slope  of  the  eddy  viscosity  is 
discontinuous  at  the  matching  point  which  separates  the  inner  and  outer 
regions,  the  Gauss-Legendre  quadrature  may  not  be  suitable  for  direct 
evaluation  of  the  integral  over  the  element  which  contains  the  matching 
point.  When  the  local  Reynolds  number  is  not  high,  the  term  involving 
the  discontinuous  slope  of  the  effective  eddy  viscosity  is  much  smaller 
than  other  terms  on  the  right  hand  side  in  equation  (3-30).  The  errors 
involving  the  integral  in  using  the  Gauss-Legendre  quadrature  on  the 
discontinuous  function  may  not  be  significant.  However,  when  the  local 
Reynolds  number  reaches  a certain  number,  the  magnitude  of  the  dis- 
continuous function  is  pronounced.  The  errors  introduced  by  this 
inaccurate  approximation  of  the  integral  involving  the  discontinuous 
function  result  in  a very  small  integration  step  in  which  the  assigned 
computer  time  is  used  up.  The  termination  of  the  computation  in  case  VI 
resulted  from  this  improper  integration.  Hence,  a special  measure  is 
required  for  the  evaluation  of  the  integral  involving  a discontinuous 
function.  For  instance,  the  integral  can  be  divided  into  two  parts  and 
then  Simpson's  rule  can  be  used  for  the  evaluation. 

In  the  process  of  computation  the  matching  point  has  to  be  deter- 
mined at  each  streamwise  station.  The  following  scheme  was  used  to  find 
the  matching  point.  First  each  element  is  further  divided  into  10  equal 
subintervals.  The  effective  eddy  viscosity  is  computed  from  the  wall 
outward  at  each  joint  node  of  the  subintervals  and  compared  with  that 
in  the  outer  region.  If  at  a certain  joint  node  the  effective  eddy  vis- 
cosity is  larger  than  that  in  the  outer  region,  the  effective  eddy 
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viscosity  is  then  computed  at  a location  which  is  less  than  the  previous 
location  by  half  the  previous  increment.  If  the  newly  computed  effective 
eddy  viscosity  is  less  than  that  in  the  outer  region,  the  effective 
eddy  viscosity  is  again  computed  at  a new  location  which  is  larger  than 
the  previous  location  by  half  the  increment.  This  procedure  is  repeated 
until  enough  points  in  the  inner  region  are  generated.  Finally,  the 
matching  point  is  obtained  by  extrapolation  using  a polynomial  of  degree 
two.  After  the  matching  point  is  determined,  Simpson's  rule  can  be  ap- 
plied to  evaluate  the  integral  over  the  special  element  in  which  the 
integrand,  is  a discontinuous  function. 

The  next  numerical  experiment  was  to  investigate  the  accuracy  of 
the  algorithm  for  determining  the  matching  point.  The  finite  element 
model  employed  to  start  the  initial  value  problem,  equation  (5-8)  and 
equation  (7-2),  is 

case  VIII:  H = 3.5,  N = 11,  h.  = /2*0.05,  0.1,  0.2,  2*0.3,  5*0.5/. 

4 

The  initial  condition  is  applied  at  Re  = 5.5  x 10  and  the  integration 

A 

routine  employed  is  the  newly  acquired  subroutine  package  LS0DE  which  is 
the  latest  version  of  the  Gear  method  developed  at  Lawrence  Livermore 
Laboratory.  The  computed  matching  point  shows  irrational  behavior  after 

c 

Re  = 1.52  x 10  . This  results  in  a small  integration  step  and  the  com- 

A 

putation  hardly  proceeds  downstream.  The  typical  jumping  characteris- 
tics of  the  computed  matching  points  are  given  in  Figure  7.  Since  the 
matching  point  represents  the  growth  of  the  inner  region,  it  should  be 
monotonically  increasing  as  the  computation  proceeds  downstream.  The 
incorrectly  determined  matching  point  may  be  due  to  improper  extrapola- 
tion or  other  causes. 
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If  the  incorrectly  determined  matching  point  is  due  to  improper 
extrapolation,  the  difficulty  can  be  overcome  by  modifying  the  original 
algorithm  for  finding  the  matching  point.  The  modified  algorithm  is 
that  one  of  the  extrapolation  points  would  be  close  enough  to  the 
effective  eddy  viscosity  in  the  outer  region  by  the  difference  of  ).5 
when  the  effective  eddy  viscosity  in  the  outer  region  is  greater  than 
50.  If  the  computed  eddy  viscosity  is  less  than  50  the  original  scheme 
is  used. 

Another  numerical  experiment  was  conducted  and  the  result  obtained 
indicated  that  the  alternate  jumping  characteristics  of  the  determined 
matching  points  is  not  due  to  inaccurate  extrapolation  but  due  to  the 
wiggle  in  the  computed  effective  eddy  viscosity.  The  wiggle  can  result 
in  difficulty  in  finding  the  matching  point  at  high  local  Reynolds 
number,  which  is  important  to  the  solution  method. 
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Figure  1.  Computed  local  skin-friction  coefficients  for  turbulent  flow  over  a flat  plate. 
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Figure  2.  Computed  effictive  eddy  viscosity  distribution 
from  case  I at  different  local  Reynolds  numbers, 
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Figure  3.  Computed  effective  eddy  viscosity  distribution  from  case  III  at 
different  local  Reynolds  numbers. 
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Figure  4.  Computed  local  skin-friction  coefficients  from  case  III  and  case  IV. 
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Figure  5.  First  derivative  of  solution  profile  with  respect  to  n at 
n=0. 05  from  case  VI . 
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Figure  6.  Second  derivative  of  solution  profile  with  respect  to 
at  n=0.0  from  case  VI. 
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Figure  7.  Jumping  of  matching  point  in  the  vicinity  of  Re  =1 . 52x1 06 

from  case  VIII: ,Rex=l .52193xl06;  — , Re  = 

1 . 52196x10  . X 


CHAPTER  VIII 
RESULTS  AND  DISCUSSION 


The  preliminary  numerical  experiments  show  that  the  choice  of  the 
first  element  size  h-j  and  other  element  sizes  in  the  inner  region  can 
be  decisive  in  the  computation  of  turbulent  flows.  The  improper  discre- 
tization of  the  flow  region  can  result  in  inaccurate  solution  profiles 
which  deteriorate  the  accuracy  of  the  effective  eddy  viscosity.  The 
wiggle  occurring  in  the  effective  eddy  viscosity  distribution  across  the 
boundary  layer  can  result  in  difficulties  in  finding  the  matching  point, 
which  is  very  important  to  the  solution  method.  Moreover,  the  properly 
chosen  value  of  H can  be  crucial  to  the  accuracy  of  the  solution  as 
well  as  to  the  efficiency  of  the  method.  The  stiffness  problem  of 
the  reduced  equations  has  been  overcome  by  using  Gear's  method. 

According  to  the  findings  in  case  VI,  the  small  first  element  size 
less  than  0.05  seems  to  be  needed  for  computations  which  proceed  to 
high  local  Reynolds  number.  The  following  two  finite  element  models 
are  designed  to  find  how  the  solution  is  affected  by  the  small  first 
element  size  as  well  as  the  element  sizes  in  the  inner  region  at  high 
Reynolds  numbers: 

case  IX:  H = 6.5,  N = 18,  h.  = /2*0.02,  2*0.03,  2*0.05,  0.2,  7*0.3, 

2*0.5,  1 .0,  2.0/ 

case  X:  H = 6.5,  N = 16,  h.  = /2*0.02,  2*0.03,  2*0.05,  0.2,  2*0.3, 

5*0.5,  1.0,  2.0/. 
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The  first  element  size  used  in  both  cases  is  0.02.  The  difference 
between  these  two  models  is  only  in  the  element  sizes  near  the  matching 
point  but  still  in  the  inner  region.  The  integration  of  the  initial 
value  problem  for  case  IX  and  case  X is  started  at  Re  = 2.2  x 106  and 

A 

c 

1.1  x 10  , respectively.  The  computed  local  skin-friction  coefficients 
in  Figure  8 compare  satisfactorily  with  the  measured  data.  The  computed 
effective  eddy  viscosity  is  therefore  found  to  be  correct.  No  wiggle 
occurs  in  the  effective  eddy  viscosity  distribution  in  either  case  in 
the  course  of  integration  up  to  Re  = 1.0  x 10  . This  indicates  that 

A 

the  wiggle  which  deteriorates  the  effective  eddy  viscosity  distribution 
can  be  due  to  the  large  first  element  size.  The  relatively  large  element 
sizes  away  from  the  boundary  may  not  be  critical  to  the  solution. 

For  the  high  local  Reynolds  number,  the  solution  profiles  change 
very  rapidly  in  the  vicinity  of  the  wall.  The  small  first  element  size 
is  absolutely  necessary  to  take  care  the  variation  of  the  solution  pro- 
file there.  So  far  we  know  that  the  first  element  size  0.02  gives  good 
results  up  to  Re  = 1 .0  x 10^,  we  would  like  to  know  whether  this  first 
element  size  0.02  is  still  valid  further  downstream.  The  finite  element 
model  for  this  purpose  was  given  by 

case  XI:  H = 4.0,  N = 15,  h.  = /2*0.02,  2*0.03,  2*0.05,  0.2,  2*0.3, 

6*0.5/. 

The  initial  condition,  equation  (7-2),  is  employed  to  start  the  initial 

4 

value  problem  at  Re  = 5.47  x 10  . The  growth  of  the  boundary  layer 

A 

thickness  is  automatically  adjusted  by  the  empirical  formula.  Whenever 
the  boundary  layer  thickness  (H)  changes,  the  increment  is  aH  = 1.  The 
number  of  elements  is  increased  by  two;  each  has  the  element  size  0.5. 


55 


The  computed  local  skin-friction  coefficient  shown  in  Figure  9 is  in 
excellent  agreement  with  the  experimental  data  up  to  Rex  = 5.0  x 107. 

The  deterioration  of  the  effective  eddy  viscosity  again  takes  place 
somewhere  after  Re  = 2.0  x 107.  The  computation  is  terminated  at 

A 

Re  = 8.6  x 107  with  H = 12  and  N = 31.  The  CPU  time  spent  on  this  case 

X 

is  250  seconds.  The  matching  point  as  a function  of  Rex,  given  in 
Figure  10,  shows  the  growth  of  the  inner  layer.  The  matching  point 
after  3 x 107  starts  to  decrease  as  the  local  Reynolds  number  increases. 
These  results  reflect  that  the  solution  profile  can  be  affected  by  the 
wiggle  in  the  effective  eddy  viscosity  after  Re  = 3 x 107.  The  computed 
displacement  boundary  layer  thickness  grows  very  fast  as  the  local  Rey- 
nolds number  increases.  Consequently,  a special  measure  in  choosing  the 
value  of  H is  needed  for  sufficiently  high  local  Reynolds  number.  The 
curvature  Wg'  and  the  first  derivative  Wg  are  also  given  in  Figures  12 
and  13,  respectively.  Both  Wg'  and  w-j  deteriorate  near  the  trailing 
edge  of  the  curves.  This  suggests  that  a very  small  first  element  size, 
less  than  0.02,  must  be  used  further  downstream. 

Since  the  accuracy  of  the  solution  profile  is  mostly  determined 
by  the  first  element  size,  the  first  element  size  obtained  by  cutting 
h.  = 0.02  of  case  XI  in  half  is  used  in  the  following  finite  element 
model  to  investigate  how  the  accuracy  of  the  solution  improved: 

case  XII:  H = 4.0,  N = 16,  h.  = /2*0.01,  0.02,  2*0.03,  2*0.05,  0.2, 

2*0.3,  6*0.5/. 

3 

The  initial  value  problem  is  solved  in  the  region  8.1  x 10  < Rex 

£ 7.8  x 107.  The  computed  skin-friction  coefficient  in  Figure  14  is 
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indeed  improved  and  is  good  up  to  7.8  x 107.  It  is  clear  that  the  small 
first  element  size  is  certainly  needed  for  accurate  prediction  of  the 
local  skin-friction  coefficient.  The  computed  effective  eddy  viscosity 
without  using  transverse  intermittency  at  several  local  Reynolds 
numbers  are  given  in  Figures  15,  16,  17  and  18,  respectively.  The  dotted 
lines  in  Figures  15,  16,  17  and  18  represent  the  effective  eddy  viscosity 
computed  by  the  inner  region  formulation  above  the  matching  point.  The 
solid  lines  represent  the  effective  eddy  viscosity  distribution  across 
the  boundary  layer.  The  computed  effective  viscosity  in  Figure  15,  16 
and  17  are  good  but  not  so  in  Figure  18.  The  wiggle  occurring  in  the 
inner  region  may  not  deteriorate  the  solution  until  the  wiggle  is  of 
sufficient  magnitude  to  cause  the  sign  of  the  first  derivative  of  the 
effective  eddy  viscosity  with  respect  to  n to  alternate.  In  order  to 
eliminate  the  wiggle  a first  element  size  smaller  than  0.01  seems  to 
be  required  further  downstream.  The  final  boundary  layer  thickness  is 
H = 12  with  24  elements.  The  increment  of  the  boundary  layer  thickness 
is  AH  = 2.  The  corresponding  increment  of  the  elements  is  2;  each  has 
element  size  1.0.  The  total  CPU  time  for  this  case  is  350  seconds. 

Because  the  computation  was  carried  out  in  the  transformed  plane, 
one  would  like  to  know  what  the  corresponding  physical  plane  looks  like. 
The  relationship  between  n and  y^  is  of  major  interest.  The  functional 
relation  obtained  from  case  XI  is  given  in  Figure  19.  The  transformation 
used  to  stretch  the  coordinate  y^  seems  to  be  suitable  for  the  boundary 
layer  flow  problems  in  this  study.  Consequently,  the  number  of  nodal 
points  required  for  the  solution  method  can  be  reduced  without  loss  of 
accuracy.  This  leads  to  computational  efficiency. 
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Figure  8.  Computed  local  skin-friction  coefficients  from  case  IX  and  case 
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Figure  9.  Computed  local  skin-friction  coefficient  from  case  XI. 


4.0 


59 


C 


Figure  10.  Growth  of  matching  point  from  case  XI. 
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Figure  11.  Growth  of  displacement  boundary  layer  thickness,  equation  (2-12), 
from  case  XI . 
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Figure  12.  Second  derivative  of  solution  profile  with  respect  to 
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Figure  13.  First  derivative  of  solution  profile  with  respect  to 
n=0.02  from  case  XI. 
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Figure  14.  Computed  local  skin-friction  coefficient  from  case  XII. 
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Computed  effective  eddy  viscosity  distribution  at 
Re  =1x1 0^  from  case  XII. 
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Figure  16. 
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Figure  17.  ComputecLeffecti ve  eddy  viscosity  distribution  at 

Re  =lxlO/  from  case  XII. 
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Figure  19.  Relationship 


CHAPTER  IX 
CONCLUDING  REMARKS 


In  this  study  the  applicability  as  well  as  the  accuracy  and 
efficiency  of  a finite  element-differential  method  has  been  investigated 
in  detail  for  steady,  two  dimensional,  incompressible  turbulent 
boundary  layer  flow  problems.  The  two-layer  eddy  viscosity  closure 
model  Wcls  utilized  to  simulate  the  Reynolds  stress  and  provide  a reason- 
able modeling  of  turbulent  flows.  The  resulting  system  of  partial 
differential  equations  has  been  transformed  into  the  proper  form  for 
application  of  the  finite  element-differential  method. 

In  the  solution  method,  the  transformed  partial  differential 
equation  is  first  reduced  into  a system  of  first  order  nonlinear 
ordinary  differential  equations  by  a subdomain  collocation  method,  in 
which  the  unknown  function  at  a streamwise  station  is  represented  by 
a classical  spline  function.  The  reduced  initial  value  problemwas  then 
integrated  numerically  by  the  modified  Hamming's  4th  order  predictor- 
corrector  method  as  well  as  by  the  Gear  method  for  stiff  equations. 

A number  of  numerical  experiments  have  been  conducted  on  the 
boundary  layer  problem  of  flow  over  a flat  plate  which  includes  laminar, 
transitional,  and  turbulent  flow  regions.  The  numerical  results  show 
that  the  solution  method  can  provide  highly  accurate  results  for  the 
complete  boundary  layer  flow  field.  In  addition,  the  application  of 
this  method  is  simple  and  straightforward.  For  example,  in  Figure  14 
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the  computed  local  skin-friction  coefficient  is  in  excellent  agreement 
with  the  experimental  data  for  the  local  Reynolds  number  up  to  7.8  x 107. 

The  efficiency  of  numerical  methods  for  boundary  layer  problems,  in 
general,  depends  on  how  many  nodal  points  are  needed  in  the  direction 
normal  to  the  flow  direction.  Due  to  the  better  approximation  of  the 
solution  profiles  by  the  spline  function  and  the  transformations  used  to 
stretch  the  y coordinate,  the  number  of  nodal  points  required  for  the 
accurate  solution  in  the  present  method  has  been  significantly  reduced  as 
compared  to  previous  methods.  For  example,  for  local  Reynolds  numbers 
less  than  2 x 10^,  the  present  method  needs  only  8 nodal  points  to  pro- 
vide highly  accurate  results,  while  the  higher  order  collocation  method 
(Rubin  and  Khosla,  1977)  and  Keller's  box  scheme  (Keller  and  Cebeci , 

1972)  may  require  20  and  30  nodal  points,  respectively.  Therefore  the 
efficiency  of  the  method  of  solution  can  be  comparable  to  finite  dif- 
ference methods. 

This  study  shows  that  when  the  local  Reynolds  number  is  sufficiently 
high,  the  solution  profiles  change  too  much  near  the  wall.  Consequently, 
due  to  this  rapid  changing  of  flow  characteristics  the  choice  of  first 
element  size  is  crucial  to  the  accuracy  of  the  solution.  The  valuable 
information  about  the  variations  of  the  first  derivative  w^  at  n-| , the 
second  derivative  Wg'  at  n0>  and  the  displacement  boundary  layer  thick- 
ness as  a function  of  the  local  Reynolds  number  has  been  obtained.  This 
information  can  give  us  a guide  to  properly  choose  the  first  element 
size.  The  results  of  numerical  experiments  show  that  a first  element 
size  of  h-|  = 0.1  is  small  enough  for  the  range  of  Wg'  less  than  1.33; 
however,  a much  smaller  value  of  h^  = 0.05  is  required  for  Wg'  between 
1.33  and  9.  For  values  of  Wg'  greater  than  9,  it  has  been  found 
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necessary  to  further  cut  down  the  first  element  size.  For  instance  the 
values  of  h]  = 0.01  and  h-j  = 0.05  will  be  needed,  respectively,  for 
24  < w“  < 32  (3.0  x 107  < Rex  < 7.0  x 107)  and  32  < w“  < 62  (7.0  x 
107  <.  Rex  <.  5.0  x 108).  For  Wq  greater  than  62,  it  is  expected  that 
the  first  element  size  which  is  less  than  0.005  should  be  used. 

The  results  obtained  in  this  study  show  that  the  finite  element- 
differential  method  can  be  efficient  and  provide  highly  accurate  results 
for  the  turbulent  flow  problems.  Moreover,  the  computer  program  can  be 
implemented  to  automatically  select  the  value  of  the  first  element  size 
based  on  the  computed  value  of  w^'  and  w-j . It  is  expected  that  there 
will  be  little  difficulty  in  extending  the  method  of  solution  to  flow 
past  more  complex  geometries  with  the  two-layer  eddy  viscosity  model 
or  with  other  transport-equation  closure  models. 


APPENDIX  A 
EDDY  VISCOSITY  MODEL 


Over  the  years,  the  eddy-viscosity  model  has  been  one  of  the  most 
popular  and  extensively  used  models.  Boussinesq  (1877)  was  the  first 
to  attack  the  problem  of  finding  a model  for  the  Reynolds  stress  -pu'v' 
by  introducing  the  concept  of  eddy  viscosity.  He  assumed  that  turbulent 
stresses  act  like  viscous  stresses,  which  implies  that  turbulent  stresses 
are  proportional  to  the  mean  velocity  gradient.  The  coefficient  of 
proportionality  was  called  "eddy  viscosity"  and  was  defined  by 


When  the  model  is  applied  to  boundary  layer  flow  problems,  the  boundary 
layer  is  regarded  as  a composite  layer  consisting  of  inner  and  outer 
regions,  and  the  distributions  of  kinematic  eddy  viscosity  are  described 
by  a separate  empirical  expression  in  each  region. 

The  determination  of  eddy  viscosity  in  the  inner  region  depends  on 
the  Prandtl's  mixing  length  theory.  The  distribution  of  the  eddy 
viscosity- in  this  region  is  given  by  (e.g.  Schlichting,  1979) 


v 


_ -pu  v 
t 3u 

p5y 


(A-l ) 


(A-2) 


where  t is  the  mixing  length  and  is  proportional  to  y,  that  is, 
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in  which  k is  the  so-called  universal  constant  and  its  numerical  value 
is  to  be  determined  from  the  experiment.  For  flows  at  high  Reynolds 
number,  k has  a value  of  0.4.  However,  in  the  viscous  sublayer,  which 
is  very  close  to  the  wall,  the  mixing  length  theory  is  not  valid.  There 
have  been  numerous  attempts  to  extend  equation  (A-3)  to  include  the 
viscous  sublayer  by  multiplying  the  mixing  length  by  some  functions. 

Van  Driest  (1956)  deduced  from  an  analysis  based  upon  experimental  data 
and  the  second  problem  of  Stokes  (1901)  that  the  correct  form  for  the 
mixing  length  in  the  inner  layer  including  the  viscous  sublayer  is 


l = ky  [1  - exp  ( - J)]  , (A-4) 

where  the  exponential  term  is  due  to  the  damping  effect  of  the  wall  on 
the  turbulent  fluctuations  and  approaches  zero  at  the  outer  edge  of  the 
viscous  sublayer  so  that  the  law  of  the  wall  as  expressed  in  equation 
(A-3)  is  valid. 

The  damping  constant  A in  equation  (A-4)  is  given  by 


A = A+  v/u  , (A-5) 

T 

where  A+  denotes  an  empirical  constant  equal  to  26  for  flat-plate  flow 
(Van  Driest,  1956)  and  u , the  friction  velocity,  was  defined  by 

T 

UT  H (xw(x)/p)i  . (A-6) 

For  incompressible  flow  with  pressure  gradient  A+  is  expected  to  vary 
somewhat  with  pressure  gradient  and  was  given  by  (Reynolds  and  Cebeci , 
1978) 


A+  = 26(1 


v dp^ 
dx 

T 


(A-7) 
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Equations  (A-4) , (A-5)  and  (A-7)  incorporated  with  equation  (A-2)  give 
the  formulation  of  the  kinematic  eddy  viscosity  distribution  in  the  in- 
ner region. 

The  formulation  of  the  eddy  viscosity  in  the  outer  region  is  based 
upon  the  Clauser  (1956)  model.  The  eddy  viscosity  in  the  outer  region 
is  scaled  by  the  displacement  boundary  layer  thickness  6*  and  the 
external  velocity  U(x).  Therefore, 

(vj  = k~U(x)<5*  . (A-8) 

z 0 c 

The  value  of  k2  in  equation  (A-8)  must  be  deduced  from  experiment  and 
has  a value  of  0.0168  (Matsui , 1963).  However,  Mellor  and  Gibson  (1966) 
used  0.016  in  their  computation  rather  than  0.0168. 

To  improve  the  calculation  of  the  shear-stress  distribution  across 
the  boundary  layer,  equation  (A-8)  needs  to  be  modified  by  a transverse 
intermittency  factor  y,  that  is, 

(vj  = k9  U (x ) <5*y  , ( A— 9 ) 

c o ^ 

in  which  y was  obtained  by  curve  fit  of  measured  data  and  given  by 
(Reynolds  and  Cebeci , 1978) 

y = [1  + 5.5  (|)6]_1  , (A- 10) 

where  6 is  the  boundary  layer  thickness. 

Equations  (A-2)  and  (A-9)  may  be  multiplied  by  a transitional  inter- 
mittency factor  which  accounts  for  the  transitional  region  between 
the  laminar  and  turbulent  flow  regions.  According  to  Chen  and  Thyson 
(1971),  the  transitional  intermittency  factor  yt  derived  from 


75 


experimental  data  for  incompressible  flow  is 


Yt 


= 1 


exp[- 


U3(x) 

1200v2 


-1.34 


(A-l 1 ) 


U(x)X 

where  x.  is  the  location  of  the  start  of  transition  and  R = . 

tr  xtr  v 

According  to  the  eddy  viscosity  model,  the  two  separate  expressions 

for  the  distribution  of  kinematic  eddy  viscosity,  including  transverse 

and  transitional  intermittency  factors  become 


(vt) _ = {ky  [1  - exp  (-  J)]  } 


9U 


(A-l 2) 


(v.)  = k7U(x)6*YYt  (A-l 3 ) 

z 0 

Equations  (A-l 2 ) and  (A- 1 3 ) can  be  used  to  describe  the  distribution  of 
the  kinematic  eddy  viscosity  for  boundary  layer  flow  with  and  without 
a pressure  gradient.  The  location  of  the  matching  point  separating  the 
inner  and  outer  regions  can  be  determined  with  the  condition 


(vt)  = (v.)  . (A-l 4 ) 

i 1 0 


The  eddy  viscosity  model  has  been  used  with  considerable  success  to 
compute  a wide  range  of  turbulent  boundary  layer  flows.  This  model 
does  not,  however,  allow  for  the  transport  of  turbulent  properties 
which  do  not  appear  parametrically  in  the  equations  of  mean  motion  and 
its  applicability  is  limited  to  near-equil ibrium  flows,  which  encompass 
most  engineering  problems. 
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APPENDIX  B 


**********  AITKEN  INTERPOLATING  FUNCTION  ********** 

FUNCTION  ATKN(  X.  V.N.K,  XI) 

AT  KN  AITKEN  INTERPOLATING  FUNCTION 

**********  AITKEN  INTERPOLATING  FUNCTION  ********** 

X - TABLE  OF  INDEPENDENT  VARIABLE  VALUES, 

(MAY  BE  ASCENDING  OR  DESCENDING), 

Y - TABLE  OF  DEPENDENT  VARIABLE  VALUES, 

(N=NO.OF  POINTS  IN  TABLES  X AND  Y.) 

K - DEGREE  CF  INTERPOLATION  DESIRED. 

XI-  X-VALUE  FOR  WHICH  INTERPOLATION  IS  DESIRED, 

THE  INTEPOLATED  VALUE  IS  RETURNED  AS  THE  FUNCTION 
VALUE 

31  CELLS  OF  BLANK  COMMON  ARE  USED. 

IMPLICIT  REAL* 8 ( A— H ,0— Z ) 

DIMENSION  X(N),Y(N),  XX  ( 13  ) • YY  ( 13) 

DATA  KMAX/  12/ 

IF  ( K ,GT.  KMAX  ,0R.  K.LE.  0 ) GO  TO  300 

K1=K+ 1 

IF(X(N)  — X(  1))  100,10,10 

10  IF  (XI-X(l))  20,20,30 
20  LL=0 

GO  TO  200 

30  IF  (X(N)-XI)  40,40,50 
40  LL=N-K1 

GO  TO  200 
50  LL— 1 
LU=N 

60  I F (LU—LL  — 1 I 180,180,70 
70  LI  = (LL+LU>/2 

IF  (X(LI)-XI)  80,80,90 
8 0 LL=H 

GO  TO  60 
90  LU=LI 

GO  TO  60 

100  IF  (XI-X(D)  120,20,20 
120  IF  (X(N)-XI)  130,40,40 

130  LL  = 1 
LU-N 

140  IF  ( LU—LL—  1 ) 180,180,150 
150  LI  = (LL*LU)  /2 

IF  (X(LI)-XI)  160,  170,  170 
160  LU=LI 

GO  TO  140 
170  LL  =L 1 

GO  TO  140 
180  LL=LL-(K 1*11/2 

IF  (LL)  20,200,190 
190  IF  (LL+K1-N)  200,200,40 
200  DO  210  1=1, K1 

I1=LL+I 

XX<  I ) = X(  I l )— XI 
210  YY ( I )=Y(  II  ) 

DO  220  I =1 , K 
DO  220  J=I,K 

22  0 YY( J*1 > = < 1 ,000/ ( XX ( J + l )— XX(  I)  ) )*( YY ( I ) *XX( J+ 1 ) — YY(  J+l  ) 
l*XX(  I ) ) 

ATKN=  YY ( K 1 ) 

RETURN 

300  *RITE<6, 30 1 ) K 

301  FORMA  T ( • 0 • , • PCLY  NO  M I AL  OF  DEGREE  K =*,I3,  2X,  'IS 

1INC0RRECT  FOR  THE  FUNCTION  SUBPROGRAM  A TK  N • ) 

RETURN 

END 
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